%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Marginal Effect for model 6: with unobs and school dummies %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%% To study the marginal effect, consider an individual in a specific
%%%% group. So first, find a group
D=100;

bnew1=bnew;
gg=71;  %% Find the groups whose mean characteristics is closest to the population.

unobs=dat(gg).unobs;

X_f = dat(gg).xt;

Mg = dat2(gg).Mg;

n=size(X_f,1);
mr=n;

k=17;

bnew=bnew1;
W_f= dat(gg).wt;
wt=W_f;
    indp_b1f=zeros(mr,1); 
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 

    for j=1:D
    
        bbf=X_f*bnew(1:k+126)+wt*X_f(:,1:k)*bnew(k+126+1:2*k+126)+Mg(:,j)*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.

        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
     
    end;
    temp3=indp_b1f/D;
    
    F=temp3;  % the F average over all 100 simulations.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 naivemarg=[2*diag(F.*(1-F))*(bnew(1)+.2*X_f(:,1)*bnew(17)) 2*F.*(1-F)*[bnew(2:16)' bnew(18+126:35+126)']];  
% not included age_2, constant as separate X here. But still include it in WX, but not report the results of its marginal effect (the 2nd to the last one).
naivemargm=mean(naivemarg(1:n,:))'*100;  

ini_n=zeros(n,k-1);

X_fa=X_f; 
for member=1:n  

    %%%%% Marginal effect for the discrete variables

    for kk=3:k-1   %%%to k-1 ,as the last regressor is squared-age
        Mg1=Mg;   
        X_fa(member,kk)=0;
                 
    indp_b1f=zeros(mr,1); 
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 
    
         for j=1:D
                
        bbf=X_fa*bnew(1:k+126)+wt*X_fa(:,1:k)*bnew(k+126+1:2*k+126)+Mg1(:,j)*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
       
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
      
        end;
      
    temp3=indp_b1f/D;
    
    F_before=temp3;  % the F average over all 100 simulations.

         X_fa(member,kk)=1;
         
 
           Mg1=Mg; 
  
    indp_b1f=zeros(mr,1);
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 
     
         for j=1:D
                
        bbf=X_fa*bnew(1:k+126)+wt*X_fa(:,1:k)*bnew(k+126+1:2*k+126)+Mg1(:,j)*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
       
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
        
        end;
        
    temp3=indp_b1f/D;
    
    F_after=temp3;  % the F average over all 100 simulations.

        amarg=F_after-F_before;
        ini_n(member,kk)=amarg(member);
    
        X_fa=X_f;
    end
end

ini_n2=ini_n(:,3:k-1);
ma1_n=mean(ini_n2)'*100; %By multiply by 100, express in terms of percentage point
naivemargm2=naivemargm;
naivemargm2(3:k-1)=ma1_n;  %%naive estimation of marginal effects

    indp_b1f=zeros(mr,1); 
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 

    for j=1:D
    
        bbf=X_f*bnew(1:k+126)+wt*X_f(:,1:k)*bnew(k+126+1:2*k+126)+Mg(:,j)*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.    
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
        indp_b2f=indp_b2f+ 1./(1+exp(2*bbf));
        
        %gradient
        deri=4./(exp(2*bbf)+exp(-2*bbf)+2);
        
        deri_f=deri_f+deri;
       
    end

fderiv_tanh=deri_f/D;   % the fderiv_tanh average over all 100 simulations.

ini=zeros(n,k-1);
affect=zeros(n,k-1);
X_fa=X_f; 
for member=1:n  
    %%%%% Marginal effect for the continuous variables
    kk=2;  %%% years in school
    dXdx=zeros(size(X_f));
    dXdx(member,kk)=1;
   
    %DXdx for contextual variables (not include dummied in contextual).
    dXdxW=zeros(n,k);
    dXdxW(member,kk)=1;
    
    partial=inv(eye(n)-diag(fderiv_tanh)*bnew(2*k+126+1)*W_f)*diag(fderiv_tanh)*(dXdx*bnew(1:k+126)+W_f*dXdxW*bnew(k+126+1:2*k+126)); %%% Calculate dM/db for group g
    amarg=2*diag(F.*(1-F))*(dXdx*bnew(1:k+126)+W_f*dXdxW*bnew(k+126+1:2*k+126)+bnew(2*k+126+1)*W_f*partial); %% marginal effect for ALL members
    ini(member,kk)=amarg(member);
    amarg(member)=0;
    affect(:,kk)=affect(:,kk)+amarg;

    kk=1; %%% age
    dXdx=zeros(size(X_f));
    dXdx(member,kk)=1;
    dXdx(member,kk+16)=2/10*X_f(member,1);
    
    dXdxW=zeros(n,k);
    dXdxW(member,kk)=1;
    dXdxW(member,kk+16)=2/10*X_f(member,1);

    partial=inv(eye(n)-diag(fderiv_tanh)*bnew(2*k+126+1)*W_f)*diag(fderiv_tanh)*(dXdx*bnew(1:k+126)+W_f*dXdxW*bnew(k+126+1:2*k+126)); %%% Calculate dM/db for group g
    amarg=2*diag(F.*(1-F))*(dXdx*bnew(1:k+126)+W_f*dXdxW*bnew(k+126+1:2*k+126)+bnew(2*k+126+1)*W_f*partial); %% marginal effect for ALL members
    ini(member,kk)=amarg(member);
    amarg(member)=0;
    affect(:,kk)=affect(:,kk)+amarg;

    %%%%% Marginal effect for the discrete variables

    for kk=3:k-1   %%%to k-1 ,as the last regressor is squared-age

        X_fa(member,kk)=0;
        
           Mg1=zeros(mr,D);

            for j=1:D;          
             
                x0=ones(mr,1)/2;
                euclid=1;

                while euclid>1e-7
                   xx_new=tanh(X_fa*bnew(1:k+126)+wt*X_fa(:,1:k)*bnew(k+126+1:2*k+126)+wt*x0*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1));
                   euclid=max(sum((xx_new-x0).*(xx_new-x0)));
                   x0=xx_new;
                end
                
            Mg1(:, j)=wt*xx_new;
            
            end
                 
    indp_b1f=zeros(mr,1);
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 

         for j=1:D
                
        bbf=X_fa*bnew(1:k+126)+wt*X_fa(:,1:k)*bnew(k+126+1:2*k+126)+Mg1(:,j)*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
              
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
     
        end;
     
    temp3=indp_b1f/D;
    
    F_before=temp3;  % the F average over all 100 simulations.

         X_fa(member,kk)=1;
       
           Mg1=zeros(mr,D);
      
            for j=1:D;          
             
                x0=ones(mr,1)/2;
                euclid=1;

                while euclid>1e-7
                   xx_new=tanh(X_fa*b1(1:k+126)+wt*X_fa(:,1:k)*b1(k+126+1:2*k+126)+wt*x0*b1(2*k+126+1)+b1(2*k+126+2)*ones(mr,1)+exp(b1(2*k+126+3))*unobs(:,j)*ones(mr,1));
                   euclid=max(sum((xx_new-x0).*(xx_new-x0)));
                   x0=xx_new;
                end
                
            Mg1(:, j)=wt*xx_new;
            
            end
                 
    indp_b1f=zeros(mr,1);
    indp_b2f=zeros(mr,1);
    deri_f=zeros(mr,1); 

         for j=1:D
                
        bbf=X_fa*bnew(1:k+126)+wt*X_fa(:,1:k)*bnew(k+126+1:2*k+126)+Mg1(:,j)*bnew(2*k+126+1)+bnew(2*k+126+2)*ones(mr,1)+exp(bnew(2*k+126+3))*unobs(:,j)*ones(mr,1); %unobs is 1 by D for each group.
         
        indp_b1f=indp_b1f+ 1./(1+exp(-2*bbf));
       
        end;
      
    temp3=indp_b1f/D;
    
    F_after=temp3;  % the F average over all 100 simulations.

        amarg=F_after-F_before;
        ini(member,kk)=amarg(member);
        amarg(member)=0;
        affect(:,kk)=affect(:,kk)+amarg;

        X_fa=X_f;
    end
end

ma1=mean(ini);  %%% average marginal effect of initiators
ma2=mean(affect)/n;%%% average marginal effect of followers
ma=[ma1; ma2]'*100%%% By multiply by 100, express in terms of percentage point
[naivemargm naivemargm2];
naive=naivemargm2






